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Charge  Density  Engineering:  A  Feasibility  Study 

Krishna  Rajan:  Department  of  Materials  Science  and  Engineering,  Iowa  State  University,  Ames  1A 
50011 

Mark  E.  Eberhart:  Department  of  Chemistry  &  Geochemistry  Colorado  School  of  Mines  Golden,  CO 
80401 

ABSTRACT 

The  objective  of  this  program  was  to  discover  composition-charge  density-property  relationships. 
The  motivation  for  this  activity  derived  from  density  functional  theory  (DFT),  which  posits  that 
all  ground  state  molecular  and  solid-state  properties  are  functionals  of  the  charge  density. 
Through  this  program  we  have  demonstrated  the  feasibility  of  coupling  advanced  ab  initio 
techniques  with  physical  insight  and  infonnatics-based  analysis  to  discover  the  nature  of  these 
relationships.  The  use  of  ab  initio  techniques  were  used  to  calculate  the  charge  density  of  metals 
and  alloys  and  reduce  the  density  to  physically  meaningful  parameters.  This  report  describes  the 
accomplishments  in  linking  ab  initio  data  with  informatics,  first  in  terms  of  charge  density  and 
secondly  in  the  related  density  of  states  (DOS).  By  uncovering  relationships  between  the 
structure  of  the  charge  density  and  the  elastic  response,  we  have  developed  chemistry-charge 
density-property  relationships.  Similarly,  we  have  (i)  identified  the  features  of  the  DOS  which 
govern  ground  state  crystal  structure  and  (ii)  “soft  modeled”  DOS  for  new  alloy  systems  without 
requiring  additional  DFT  calculations.  We  have  shown  that  the  charge  density  descriptors,  when 
analyzed  within  a  statistical  framework,  uncover  previously  unknown  relationships  between 
crystallographic  structure  and  charge  density. 


BACKGROUND 

As  is  known  from  density  functional  theory,  the  ground  state  properties  of  an  atomic  system  are 
detennined  solely  by  its  charge  density.  As  mechanical  response  at  typical  strain  rates  is 
dominated  by  ground  state  properties,  it  may  be  possible  to  relate  features  of  the  charge  density 
to  specific  mechanical  properties.  If  these  relationships  could  be  combined  with  an 
understanding  as  to  the  effect  of  elemental  composition  on  these  electronic  features,  one  could 
manipulate  charge  density  so  as  to  produce  desirable  properties.  That  is,  one  could  design  the 
chemistry  of  solids  to  yield  optimum  mechanical  response.  This  program  was  designed  to  test 
these  speculations.  Specifically  we  explored  two  sets  of  relationships:  first,  a  set  that  relates  the 
geometry  of  the  charge  density  to  a  material’s  elastic  constants  (a  measure  of  a  purely  ground 
state  mechanical  property),  and  second,  a  set  that  relates  elemental  composition  to  the  structure 


of  the  charge  density.  Combined,  these  two  sets  provide  the  process-structure-property 
relationships  for  a  new  discipline  of  charge  density  engineering. 


Our  approach  builds  on  the  quantum  theory  of  atoms  in  molecules  (QTAIM),  which  has 
demonstrated  that  the  structure  of  the  charge  density  may  be  compactly  represented  by  a  few 
parameters  describing  the  density  around  specific  points.  This  “local  approximation”  reduces  the 
description  of  a  3D  scalar  field — the  charge  density — to  that  of  identifying  the  locations  of  the 
specified  points  and  computing  the  values  of  shape  descriptors  at  these  points.  This  program 
used  first  principle  calculations  to  create  databases  of  charge  density  descriptors  combined  with 
physical/chemical  arguments  to  identify  meaningful  combinations  of  these  descriptors  which  are 
then  subject  to  statistical  and  data  mining  techniques  to  establish  correlations  between  these 
meaningful  descriptors  and  measured  or  calculated  single  crystal  elastic  constants. 

In  previous  work  we  have  shown  that  the  fee  single  crystal  elastic  constant  of  simple  metals  can 
be  correlated  with  the  descriptors  around  a  single  critical  point — the  bond  point.  This  correlation 
was  demonstrated  by  inspection,  with  the  introduction  of  a  geometrically  meaningful 
representation  of  the  Hessian  tensor  in  which  the  pi  1,  p22,  and  p33  were  combined  to  yield  three 
new  tensor  invariants.  In  addition  to  this  representation  of  the  Hessian  invariants  at  a  critical 
point,  there  are  many  others,  such  as  the  Gaussian  curvature,  G  =  pnp22p33,  and  the  mean 
curvature  or  Laplacian,  L  =  pn+p22+p33.  Attempts  to  extend  this  correlation  between  shape 
descriptors  and  elastic  constants  to  other  crystal  structures  were  ultimately  frustrated  by  the  fact 
that  the  elastic  constants  were  not  dominated  by  the  shape  parameters  around  a  single  type  of  cp 
but  around  all  cps.  In  these  cases,  we  needed  to  look  for  correlations  between  elastic  constants 
and  often  more  than  20  shape  parameters,  which  could  be  represented  in  several  different 
meaningful  forms. 

We  made  significant  progress  in  both  the  ab  initio  and  infonnatics  areas  in  the  previous  reporting 
periods,  with  the  accomplishments  of  the  two  research  thrusts  fully  integrated  in  the  final  year  of 
the  program,  as  discussed  in  this  report.  The  previously  reported  accomplishments  in  ab  initio 
included:  (i)  creation  of  a  data  set  giving  the  shape  parameters  at  all  bee  and  fee  transition  metals 
at  critical  points;  (ii)  calculated  the  electronic  structure  of  160  B2  intermetallics  and  calculated 
their  single  crystal  elastic  moduli;  (iii)  located  and  calculated  the  shape  parameters  for  ~80  of 
these  compounds;  (iv)  constructed  an  empirical  model  to  account  for  the  observed  trends  in  the 
single  crystal  B2  elastic  constants;  and  (v)  created  a  database  of  30  entries  to  provide  initial 
correlations  of  the  effects  of  composition  on  shape  descriptors.  The  previously  reported  progress 
in  infonnatics  analyses  included:  (i)  performed  a  quantitative  assessment  of  correlations  between 
topological  descriptors  from  critical  points  using  data  dimensionality  reduction  methods;  (ii) 
confirmed  that  the  chosen  shape  parameters  may  serve  as  robust  modeling  parameters  for  the 
next  stage  of  the  work  by  demonstrating  that  they  serve  as  strong  classifiers  of  crystal  structure; 
and  (iii)  captured  systematics  in  the  impact  of  transition  metal  rare  earth  chemistry  on  B2  charge 
density  shape  descriptors. 


The  goals  accomplishments  in  linking  the  informatics  and  ab  initio  aspects  of  this  program 
which  are  the  focus  of  this  final  report  are  summarized  as  follows.  First,  given  the  complexity  of 
critical  point  data,  the  challenge  is  to  link  the  shape  of  the  charge  density  topology  to  chemistry 
and  property,  which  represents  a  significant  data  challenge.  To  address  this  challenge,  we  have 
followed  a  two  prong  approach,  where  first  we  predict  elastic  moduli  as  a  function  of  the  charge 
density  topology,  which  accomplishes:  (i)  accelerated  prediction  of  elastic  constant  data,  and  (ii) 
reduced  the  amount  of  infonnation  needed  from  ab  initio  calculations.  Given  this  limited 
knowledge  requirement,  we  then  predict  the  critical  point  descriptors  based  only  on  elemental 
descriptors  from  the  periodic  table.  Therefore,  the  final  outcome  is  that  we  can  predict  critical 
charge  density  values  from  basic  elemental  properties,  and  from  these  critical  charge  density 
values  we  can  then  predict  the  mechanical  response,  thus  accelerating  the  selection  of  targeted 
chemistries.  Additionally,  we  have  studied  the  related  DOS  curves  for  both  fundamental  science 
and  engineering  purposes.  The  fundamental  aspect  of  this  approach  is  to  answer  the  century  old 
question  of:  why  do  metals  have  their  ground  state  crystal  structures?  The  engineering  objective 
is  to  predict  the  DOS  of  alloys  given  the  DOS  of  the  constituent  elements,  and  then  from  this 
informatics  alloy  DOS  to  extract  the  moduli.  These  integrated  ab  initio  and  informatics  thrusts 
serve  to  guide  future  calculations,  modify  the  chemical  search  space,  and  accelerate  design  of 
new  materials  for  targeted  electronic  behavior. 


Detailed  Progress  Report  -  Ab  Initio  (Lead  -  Colorado  School  of  Mines) 

(i)  Building  ab-initio  database 

In  addition  to  several  small  data  sets  for  the  pure  fee  and  bee  transition  metals,  we  completed  the 
construction  of  a  large  data  base  of  shape  descriptors  and  elastic  moduli  for  intermetallics  with 
the  B2  (CsCl)  crystal  structure.  We  chose  to  investigate  14  common  B2  intermetallics  and  128 
RM  (R  =  rare  earth  elements,  M  =  metallic  elements  from  groups  2,  8-13).  These  compounds 
were  chosen  to  give  a  broad  range  of  elastic  responses.  All  calculations  were  performed  with 
VASP  using  the  PBE  and  PW91  generalized  gradient  approximations.  We  proceeded  by 
generating  the  charge  density  on  a  25  point/angstrom  mesh,  locating  the  CPs  of  the  charge 
density,  and  generating  the  geometric  properties  of  the  charge  density  at  those  points,  using  our 
software  package  TECD.  Notably,  we  discovered  that  there  are  two  characteristic  B2  topologies, 
one  with  second  neighbor  B-B  bonds  only,  and  one  with  both  second  neighbor  A-A  bonds  and 
second  neighbor  B-B  bonds  (Figure  1).  In  both  cases  there  are  first  neighbor  A-B  bonds.  The 
single  type  of  second  neighbor  bond  topology  is  more  common  in  this  set  of  compounds.  Only 
the  13  materials  listed  in  Table  1  display  both  A-A  and  B-B  bonds. 


Fig  1:  The  two  bonding  topologies  of  the  B2  compounds  investigated.  The  A-B  and  B-B  bonding 
topology  is  shown  left,  with  turquoise  spheres  representing  both  first  and  second  neighbor  bond  points. 
The  A-B,  A-A,  and,  B-B  bonding  topology  is  shown  right,  where  the  additional  B-B  bond  points  are 
shown  as  blue  spheres. 


Table  1:  The  B2s  displaying  both  A-A  and  B-B  second  neighbor  bonding. 


AgHo 

HgTb 

AgNd 

RhYb 

AuCd 

ZnEr 

AuPr 

ZnHo 

CdLa 

ZnTb 

HgEr 

ZnYb 

HgNd 

In  addition  to  calculating  charge  density  we  also  calculated  the  elastic  moduli  of  all  the  B2 
crystals.  The  single  crystal  elastic  moduli,  Cn,  C12,  and  C44,  were  computed  from  the  curvature  of 
the  internal  energy  versus  strain  curves.  The  strains  that  resist  Cn,  Ci2,  and  C44  are  as  follows: 

(Cu)  8n  =  622  =  y,  £33=(1+Y)  2-1 

(C  L2)  £  1 1  =  £  22  =  £  33  =  Y 

(C44)  £[2  =  £  23  =  £31  =  Y 


In  all  cases  we  calculated  the  internal  energy  using  applied  strains  ranging  from  -0.025  to  0.025 
in  steps  of  0.005.  For  the  common  B2  intermetallics,  our  calculated  results  agree  with  published 
experimental  (calculated)  results  (Table  2).  The  independent  shear  constants  are  C44  and  C’  = 
l/2(Cj  i-C  i2)  and  the  strains  resisted  by  these  constants  are  depicted  in  Figure2. 


Fig  2:  Strains  with  non-vanishing  £xv  (left)  and  £xx  =  (right).  The  arrows  represent  the  strain  and  the 
red  dotted  lines  the  nodes  of  that  strain. 


Table  2:  A  comparison  of  ab  initio  elastic  properties  of  selected  B2  intermetallics  (this  study)  with  the 
available  reported  data. 


Phase 

C„(GPa) 

Ci2(GPa) 

C44(GPa) 

A 

CuZr 

139.2 

111.4 

44.7 

3.604 

137.3 

108.9 

44.9 

3.160 

138.0 

112.0 

45.0 

3.462 

TiZn 

140.7 

107.2 

97.0 

5.791 

140.6 

106.2 

99.2 

5.767 

143.5 

99.0 

94.0 

4.225 

NiTi 

176.5 

156.1 

46.8 

4.582 

179.5 

156.4 

49.5 

4.291 

183.0 

146.0 

46.0 

2.487 

178.2 

147.6 

49.0 

3.202 

AlZr 

145.4 

85.2  86.2 

27.1 

0.901 

145.7 

29.8 

1.000 

CuZn 

121.7 

110.3 

80.7 

14.211 

123.4 

110.9  116.0 

84.3 

13.480 

124.0 

79.0 

AgMg 

80.9 

56.6 

47.4 

3.893 

79.8 

55.5 

47.8 

3.930 

AgZn 

97.1 

88.7 

51.8 

12.40 

97.3 

88.0 

54.0 

11.54 

AgCd 

78.3 

74.9 

44.5 

25.77 

76.8 

72.7 

45.1 

22.27 

AuZn 

123.7 

113.9 

44.0 

8.969 

124.1 

113.5 

46.8 

8.830 

AuCd 

93.2 

91.6 

38.4 

45.17 

91.8 

89.9 

39.6 

42.75 

AlTi 

150.3 

96.2 

68.4 

2.526 

150.9 

97.6 

70.1 

2.630 

NiAl 

210.1 

135.8  116.4 

115.9 

3.127 

207.0 

106.8 

3.172 

PtTi 

203.8 

181.2 

46.9 

4.137 

203.3 

182.2 

48.5 

4.592 

PdTi 

169.8 

150.3 

44.7 

4.585 

168.6 

150.4 

45.9 

5.049 

CuY 

115.0 

47.6  47.1  47.4 

35.8 

1.062 

114.3 

47.7 

36.7 

1.065 

116.4 

34.5 

1.000 

116.0 

31.9 

0.934 

AgY 

98.1 

53.3  51.9  50.3 

35.5 

1.587 

96.1 

54.5 

35.7 

1.615 

105.3 

37.2 

1.353 

102.4 

32.6 

1.361 

Of  the  142  B2  intermetallics  for  which  the  electronic  structure  was  calculated,  we  were  able  to 
obtain  reliable  geometric  charge  density  descriptors  for  73  (Table  3).  Our  qualitative 
understanding  is  rooted  in  detennining  how  the  charge  density  responds  to  an  applied  strain,  i.e. 
the  nature  of  the  charge  flow  or  redistribution.  Ideally,  we  would  describe  this  charge  flow  using 
bond  bundles  [1,2].  The  B2  structure  is  characterized  by  only  two  types  of  bond  bundles  and  as 
charge  flows  from  one  it  must  flow  into  the  other,  which  makes  the  analysis  of  quantities  that  go 
as  the  ratio  of  bond  bundle  properties  simpler  to  rationalize  than  those  related  to  the  absolute 
value  of  bb  properties.  One  such  property  is  the  elastic  anisotropy,  A  =  2C44/(Cn-Ci2).  Hence, 
we  began  our  development  of  a  qualitative  understanding  of  elastic  response  by  attempting  to 
correlate  elastic  anisotropy  with  critical  point  shape  descriptors.  Figure  3  shows  the  bond  points 
of  the  common  B2  structure.  The  crystal  is  composed  of  8  first  neighbor  bond  bundles,  one 
connecting  the  central  atom  to  each  vertex,  and  6  second  neighbor  bond  wedges,  one  connecting 
the  central  atom  to  each  face.  By  symmetry  the  nuclei  of  the  bond  bundles  are  pinned,  along  with 
the  cage  points.  Ring  points  can,  however,  move  in  the  (100)  plane  in  the  <11 0>  direction,  as 
shown  in  the  lower  panes  of  Fig.  3.  Only  movement  of  the  ring  point  changes  the  volume  of  the 
bond  bundles.  We  can  then  use  the  distance  between  the  ring  point  and  the  first  neighbor  bond 
point  as  an  indicator  of  the  volume  of  the  first  neighbor  bond  bundle.  Finally,  we  can  get  an 
approximation  of  the  charge  in  the  bond  by  multiply  this  distance  by  the  number  of  electrons  at 
the  bond  point. 


Table  3:  The  chemical  formula  of  the  B2  intermetallics  with  reliable  charge  density  descriptors. 


AgDy 

AuCd 

CdDy 

CdYb 

HgHo 

IrYb 

RhTb 

AgEr 

AuDy 

CdEr 

CuDy 

HgLa 

PdDy 

RhY 

AgGd 

AuEr 

CdGd 

CuEr 

HgNd 

PdEr 

RhYb 

AgHo 

AuGd 

CdHo 

CuGd 

HgPr 

PdHo 

ZnCe 

AgLa 

AuHo 

CdLa 

CuHo 

HgSm 

PdYb 

ZnDy 

AgNd 

AuNd 

CdNd 

CuTb 

HgTb 

RhDy 

ZnEr 

AgPr 

AuPr 

CdPr 

CuY 

HgY 

RhEr 

ZnGd 

AgSm 

AuSm 

CdSm 

HgDy 

IrEr 

RhGd 

ZnHo 

AgTb 

AuTb 

CdTb 

HgEr 

IrHo 

RhHo 

ZnLa 

AgY 

AuY 

CdY 

HgGd 

IrY 

RhSm 

ZnSm 

ZnTb 

ZnY 

ZnYb 

Fig  3:  The  edges  of  a  first  neighbor  (left)  and  second  neighbor  bond  bundle  (right)  in  a  B2  material. 
Bond  points  are  shown  as  turquoise,  rings  as  red,  and  cages  as  yellow  spheres,  while  bond  paths  are 
shown  as  gray  cylinders  in  the  top  to  panes.  The  unit  cell  is  composed  of  8  Is'  neighbor  bond  bundles 
and  6  second  neighbor  bond  wedges  (half  a  bond  bundle).  The  lower  two  panes  show  the  cross  section 
of  the  bond  bundles  with  arrows  indicating  the  directions  the  ring  points  can  move.  The 
semitransparent  ring  points  also  illustrate  these  directions. 


We  would  expect  that  as  this  rough  measure  of  charge  increases  the  elastic  anisotropy  decreases 
because  as  the  amount  of  charge  in  the  first  neighbor  bonds  increases  so  does  the  charge  flow 
from/to  those  bonds.  Fig.  4  supports  this  assumption.  It  is  a  plot  of  the  distance  between  the  ring 
and  first  neighbor  bond  point  (in  fractional  coordinates)  multiplied  by  the  charge  at  the  bond 
point  versus  the  elastic  anisotropy.  From  the  above  correlation  we  conclude  that  the  Guassian 
curvature  at  the  second  neighbor  bond  points  (cage  points)  correlates  with  C’,  but  when  we 
attempt  to  correlate  C44  with  eignevalues  of  H(cp)  (or  invariants  of  H(cp)),  we  find  no  good 
correlations.  This  observation  is  probably  due  to  the  fact  that  the  charge  flow  to/from  the  CPs 
under  strains  with  non-vanishing  £xx  =  -£„,  occurs  along  the  principal  directions,  whereas  under  a 
strain  with  non-vanishing  £xy  the  charge  flow  is  not  along  a  principal  direction.  Such  a 
phenomenon  has,  to  our  knowledge,  not  been  reported.  The  simplest  possible  measure  of  this 
flow  would  be  given  by  the  curvature  along  the  ring-bond  ridge,  as  it  is  the  direction  of  charge 
flow. 
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Fig  4:  Plot  of  normalized  distance  (distance/lattice  constant)  between  the  ring  point  and  1st  neighbor 
bond  point  multiplied  by  the  charge  density  at  the  1st  neighbor  bond  point  versus  the  calculated  elastic 
anisotropy  (2  C4/(Cu-Ci^).  The  normalized  distance  serves  as  an  approximation  of  the  volume  of  the  1st 
neighbor  bond  bundle. 


(ii)  Generalized  Bonding  Model 

Just  as  a  two-dimensional  surface  may  be  triangulated,  a  three  dimensional  volume,  such 
as  the  charge  density,  may  be  covered  by  a  set  of  irregular  tetrahedra  sharing  faces  edges 
and  vertices.  The  resulting  system  is  called  a  simplical  complex.  Being  tetrahedral,  the 
irreducible  bundle  fonns  a  simplical  complex  over  the  charge  density.  In  this  unique 
form,  the  1 -skeleton  or  underlying  graph  of  the  simplical  complex  takes  on  special 
meaning.  The  1 -skeleton  of  a  3D  simplical  complex  is  the  union  of  all  the  tetrahedral 
edges,  which  in  this  case  is  the  set  of  all  the  system’s  1 -ridges  and  by  default  includes  the 
bond-path,  which  is  a  1 -ridge  connecting  adjacent  charge  density  maxima  through  a  (3,  - 
1)  saddle  point.  Thus,  the  underlying  graph  contains  as  a  subset  the  molecular  graph 
common  to  depictions  of  the  bonding  in  molecules  and  solids.  Our  conjecture  is  that  the 
entire  underlying  graph  (1 -skeleton),  as  a  set  of  pairwise  connections  between  charge 
density  critical  points,  provides  an  excellent  quantitative  representation  of  the  bonding 
and  elastic  properties  in  a  molecule  or  solid. 

As  a  step  in  confirming  this  conjecture,  we  modeled  inter-critical  point  spring  constants 
by  investigating  crystalline  systems.  Our  first  effort  was  to  recover  the  Cauchy  pressure 
for  the  BCC  transition  metals,  knowing  only  information  about  the  charge  density  and  its 
Hessian  at  the  non-nuclear  critical  points.  We  started  with  the  nonmagnetic  BCC 
transition  metals  because,  of  the  common  crystal  structures,  only  the  BCC  structure 


possesses  the  simplest  possible  topology  with  one  symmetry  unique  nuclear,  bond,  ring 
and  cage  CP  and  hence  (the  minimum)  six  different  1 -ridges.  We  adopt  the  notation 
kcpicp2  to  denote  the  spring  constant  of  the  link  joining  cpl  to  cp2,  where  a  cp  can  be  n,  b, 
r,  or  c  for  nuclear,  bond,  ring  and  cage  respectively.  Thus,  kbr  corresponds  to  the  spring 
constant  of  the  bond  to  ring  CP  connection. 

As  with  all  cubic  crystals,  the  elastic  response  of  this  entire  set  of  springs  is  tied  to  three 
independent  elastic  constants:  Cn,  C12,  and  C44  giving: 

-1/3 

Ci  1  —  Va  ( 0.67  knb  +6.80  km-  +2  knc  +13.96  kbr  +7.31  ktc  +14.56  krc ) 

-1/3 

C12  ~  Va  {0.67  k„b  +1.60  knr  +0  knc  +5. 76  kbr  +3.31  kbc  +1.69  kIT) 

-1/3 

C44  —  Va  {0.67  knb  +1.60  knr+0  knc  + 5.67  kbr  +1  kbc  +0.80  krc ) 


where  Va  is  the  atomic  volume.  Note  that  the  quantity  Cn~  C44,  (the  Cauchy  pressure),  is 
a  function  of  only  three  spring  constants 

-1/3 

C12  -  C44  =  Va  (3.10  khr  +2.31  kbc  +0.89  krc), 
where  all  tenns  related  to  connections  with  the  nuclear  cp  have  dropped  out. 

Systems  where  the  Cauchy  pressure  is  zero  are  indicative  of,  “atoms  interacting  through  a 
central  force  potential,”  which  includes  simple  spring  like  connections  and  accounts  for 
the  fact  that  connections  with  the  nuclear  cp  have  dropped  out.  Positive  deviations  reflect 
“metallic  bonding,”  which  necessitates  corrections  involving  “many  body  interactions  that 
result  when  an  atom  interacts  with  the  electron  gas  of  its  neighbors.”  Negative  Cauchy 
pressures  derive  from  “covalent  bonding  character”  and  require  angular  corrections.  Our 
task  reduced  to  finding  relationships  between  the  shape  of  the  charge  density  at  the 
connected  critical  points  and  the  intercritical  point  spring  constants. mGeometric  and 
physical  intuition  suggested  that  the  parameters  of  interest  would  be  given  the  charge 
density  at  a  critical  point,  p,  and  an  angle,  0,  related  to  the  eccentricity  of  the  charge 
density  about  the  critical  point,  giving  four  descriptors  for  each  intercritical  point  spring 
constant. 


(Hi)  Energetic  Materials 

A  picture  of  impact  sensitivity  based  on  the  bond  bundles  of  the  electron  charge  density 
has  been  developed,  which  includes  the  response  of  both  inter-  and  intramolecular 
bonding  interactions.  Impact  sensitive  materials  were  found  to  have  closed  intramolecular 
bond  bundles  with  a  low  electron  count  that  serves  as  a  trigger  linkage,  while  insensitive 


materials  do  not.  The  shape  and  electron  count  of  the  intramolecular  bond  bundles  was 
found  to  change  between  the  gas  phase  and  solid  state  due  to  the  fonnation  of 
intennolecular  bonds.  In  the  case  of  polynitrobenzenes  this  change  was  subtle  and  did  not 
affect  the  trigger  linkages.  However,  the  intennolecular  bonds  in  crystalline  RDX 
transform  the  C-N  trigger  linkages  found  in  the  gas  phase  to  N-N  trigger  linkages  in  the 
solid  state.  This  observation  offers  a  theoretical  justification  of  the  well-known 
experimentally  observed  differences  in  the  decomposition  of  gas  phase  and  crystalline 
RDX. 


(iv)  Critical  Point  Input  into  Informatics  Analysis 


Table  4.  A  list  of  descriptors  for  electronic  charge  density  topology  of 
intermetallic  materials.  These  values  were  calculated  for  the  B2  intermetallics 
listed  in  Table  3.  These  descriptors  were  analyzed  by  informatics  for  developing 
structure-property  relationships,  where  the  property  is  based  on  the  elastic 
moduli. 


No. 

Descriptor 

Description 

1 

kl  bb 

Eigenvalue  1  at  the  bond-bundle  critical  points 

2 

kl_bb 

Eigenvalue  2  at  the  bond-bundle  critical  points 

3 

^3_bb 

Eigenvalue  3  at  the  bond-bundle  critical  points 

4 

Pbb 

Charge  at  the  bond-bundle  critical  points 

5 

tan(0)bb 

(kz  bb  /  ^l_bb)1/_ 

6 

A_bb 

Laplacian  (=ki _bb  +  k2_bb  +  k3_bb) 

7 

kl_b 

Eigenvalue  1  at  the  bond  critical  points 

8 

^2_b 

Eigenvalue  2  at  the  bond  critical  points 

9 

hb 

Eigenvalue  3  at  the  bond  critical  points 

10 

Pb 

Charge  at  the  bond  critical  points 

11 

tan(0)b 

(Wki_b)I/2 

12 

A  b 

Laplacian  (=A,!  b  +  f  b  +  hb) 

13 

ki_c 

Eigenvalue  1  at  the  cage  critical  points 

14 

^2_c 

Eigenvalue  2  at  the  cage  critical  points 

15 

^“3_c 

Eigenvalue  3  at  the  cage  critical  points 

16 

Pc 

Charge  at  the  cage  critical  points 

17 

A_c 

Laplacian  (=ki  c  +  f  c  +  ^3c) 

18 

kl_r 

Eigenvalue  1  at  the  ring  critical  points 

19 

A.2_r 

Eigenvalue  2  at  the  ring  critical  points 

20 

^3_r 

Eigenvalue  3  at  the  ring  critical  points 

21 

Pr 

Charge  at  the  ring  critical  points 

22 

A  r 

Laplacian  (=k;  ,  +  ^21  +  ^) 

23 

fl 

Electron  density  flatness  (=pcmin/ pb™ax) 

24 

f2 

Energy  term  related  to  stacking  fault  (=(pb  -  pc)2) 

25 

f*3 

Energy  term  related  to  surface  (=(pb  +  pc)2) 
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(i)  Genetic  Programming  for  QSPRs 

Genetic  programming  is  an  approach  that  has  been  applied  for  optimizations  in  many 
engineering  fields  such  as  process  control.  We  apply  it  for  the  first  time  to  ab  initio 
derived  data  to  predict  the  elastic  constants  of  transition  metals  as  a  function  of  critical 
point  shape  descriptors.  This  work  serves  two  purposes:  (i)  QSPR  can  be  used  to  screen 
for  material  properties  and  accelerate  the  identification  of  new  candidate  materials;  and 

(ii)  the  terms  in  the  QSPR  model  provide  physical  understanding  of  material  behavior  that 
otherwise  is  too  complex  to  identify. 

We  have  used  genetic  programming  [2]  to  obtain  a  function  of  the  charge  density  shape 
descriptors  which  gives  the  spring  constants  of  the  intercritical  point  connections.  GP  is 
an  evolutionary  algorithm  based  on  the  ideas  of  biological  evolution,  designed  to  evolve 
computer  programs  that  optimize  a  user-defined  function.  The  objective  is  to  create  a 
population  of  possible  solutions  to  the  problem.  The  members  of  the  population  are 
evaluated  using  a  fitness  function  (in  this  case  the  mean  of  the  squared  prediction  errors 
for  each  of  the  five  BCC  metals)  that  measures  how  good  the  functions  are  in  predicting 
the  desired  quantity  (Cauchy  pressure)  from  the  (shape)  descriptors. 

A  genetic  programming  search  was  conducted  for  a  function  that  would  express  the  spring 
constant  between  two  critical  points  as  a  function  of  four  variables  mentioned  above. 
Certain  constraints  were  imposed  on  the  function,  such  as  maintaining  symmetry  and 
simplicity  of  the  model  to  make  physical  interpretation  possible.  The  fitness  function  was 
set  to  the  mean  of  the  squared  deviations  of  the  predicted  value  of  Cauchy  pressure  from 
the  value  calculated  through  DFT.  Two  different  function  forms  were  found  that  obeyed 
the  defined  criteria  and  were  highly  accurate: 

Function  1: 


'jcp\cp2 


0.1079(f  -flcpl)(f  -ecp2) 
PcplPcp2 


Function  2: 


k'cp  lcp2  =  18.22(0.5784  -  Pcp i)(0.5784  -  pcp2)(2.063  -  0cpl)(2.O63  -  0cp2) 

Function  1  is  optimized  for  simplicity  and  predicts  the  Cauchy  pressure  with  more  than 
97%  accuracy.  A  slightly  better  accuracy  of  prediction  (98.6%)  is  possible  with  Function 


2,  though  at  the  cost  of  the  simplicity  of  the  relation  between  the  spring  constant  and  the 
critical  point  descriptors.  The  plots  for  the  predicted  and  actual  Cauchy  pressure 
calculated  through  DFT  are  presented  in  Figure  5. 


Actual  Cauchy  Pressure  Actual  Cauchy  Pressure 


Figure  5.  Genetic  programming  prediction  of  Cauchy  pressure  versus  the  DFT 
predicted  Cauchy  pressure  for  bcc  transition  metals.  The  model  optimized  for 
simplicity  has  97.3%  accuracy,  while  the  model  optimized  for  accuracy  has  98.6% 
accuracy. 

Our  analysis  of  these  two  functional  forms  provides  insight  as  to  the  origin  of  the 
intercritical  point  spring  constant.  The  Cauchy  pressure  goes  to  zero  when  the  charge 
density  takes  on  a  shape  resulting  from  the  overlap  of  spherically  symmetric  densities. 

This  “reference  state,”  though  built  from  overlapping  spherical  charge  densities,  will  likely 
be  characterized  by  a  BCC  topology  with  nonzero  values  for  all  the  shape  descriptors  of 
the  observed  systems. 

If  we  assume  the  tails  of  these  overlapping  spherical  charge  densities  go  as  e  ar.  It  is  then 
an  easy  matter  to  calculate  the  values  of  the  reference  quantities  as  functions  of  a.  While 
the  reference  quantities,  p  and  9,  vary  with  a,  kcpicp2,  as  calculated  be  either  Function  1  or 
Function  2,  vary  more  slowly  and  are  close  to  zero.  The  compelling  conclusion  is  that  the 
inter-critical  point  spring  constants  provide  information  about  the  amount  of  charge  that 
flows  from  a  CP  in  its  reference  state  to  yield  the  same  CP  in  the  observed  state;  a  number 
of  functional  forms  may  scale  with  the  “distance”  between  these  two  states. 


If  we  are  to  generalize  this  work  to  different  topologies,  we  must  identify  a  true  invariant 
of  overlapping  spherical  charge  densities,  or,  a  best  method  to  estimate  the  values  of  the 
shape  descriptors  for  the  reference  state  of  given  element  in  a  given  crystal  structure. 

(ii)  Design  Rules  for  New  Intermetallics 

The  predictions  of  mechanical  property  as  a  function  of  charge  density  has  been  extended 
to  intermetallic  compounds,  with  elastic  constants  predicted  as  a  function  of  25  shape 
descriptors.  Additionally  these  descriptors  are  used  to  predict  the  stability  of  a 
compound;  that  is,  to  predict  whether  a  new  material  could  fonn  in  a  particular  crystal 
structure.  These  results  have  been  tested  on  known  systems,  while  future  work  will 
expand  this  work  to  suggest  new  candidate  materials. 

We  found  that  the  predictive  models  achieved  by  correlating  topological  charge  density  in 
tenns  of  critical  points  with  elastic  constants  discloses  the  hidden  structure-property 
relationships  that  cannot  be  identified  from  crystallographic  information.  The  results  for 
prediction  of  the  elastic  constants  from  the  charge  density  critical  points  are  provided  in 
Figure  6.  Our  results  demonstrate  that  charge  density  topology  data  of  a  material  system 
contains  the  infonnation  required  for  the  prediction  of  elastic  constants  of  a  new  material. 
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Figure  6.  Prediction  of  elastic  constants  based  on  a  partial  least  squares 
regression  model  of  73  intermetallic  compounds.  The  accuracy  of  this  approach 
demonstrates  the  ability  to  predict  mechanical  properties  for  alloys. 


In  addition  to  predicting  mechanical  properties,  we  have  developed  a  series  of  design 
rules  for  predicting  mechanical  stability  of  new  compounds,  building  on  the  ability  to 
define  mechanical  properties  as  a  function  of  critical  point  descriptors.  For  a  cubic 
crystal,  the  following  three  conditions  are  generally  accepted  as  the  elastic  stability 
criteria: 

(i)  Bulk  modulus,  K  =  l/3(Cn  +  2C 12)  >  0 

(ii)  Shear  modulus,  G  =  C44  >  0 

(hi)  Cprime  =  l/2(Cn  -  C12)  >  0 


If  all  of  these  three  conditions  are  not  met  simultaneously,  the  cubic  crystal  is 
(mechanically)  unstable.  The  charge  density  data  were  correlated  with  the  stability 
condition  using  a  recursive  partitioning  algorithm.  As  a  result,  if-then  rules  of  stability 
criteria  were  developed  as  shown  in  Figure  7.  Only  five  critical  point  descriptors  were 
selected  for  building  the  tree  model,  indicating  that  those  five  descriptors  are  the  defining 
factors  for  the  detennination  of  stability  condition  for  B2  compounds.  Starting  from  the 
root  node  of  the  if-then  tree  structure,  the  stability  of  a  new  material  can  be  predicted. 
According  to  the  internal  validation  of  the  model,  R~=0.8526  was  obtained.  This  shows 
that  the  elastic  stability  of  a  material  can  be  predicted  directly  from  the  electronic  charge 
density  infonnation  of  the  material. 
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Figure  7.  Classification  tree  model  for  the  prediction  of  the  stability  of  B2  compounds  in 
terms  of  the  critical  points  of  topological  charge  density;  the  numbers  shown  at  the  end  of 
each  branch  indicate  the  number  of  potential  new  compounds  classified  into  each  branch. 


(in)  Identification  of  Minimal  Electronic  Structure  Information 


One  of  the  primary  objectives  of  informatics  is  to  identify  the  minimal  amount  of  information 
necessary  for  describing  a  material.  To  this  end,  we  have  applied  principal  component  analysis 
(PCA)  and  variable  importance  parameter  (VIP)  in  identifying  the  parameters  which  are  most 
important  in  describing  electronic  structure.  By  identifying  the  fewest  parameters  needed  to 
model  the  elastic  moduli  from  the  charge  density  topology,  the  infonnatics  modeling  of  these 
parameters  becomes  manageable,  as  discussed  in  the  next  section. 

To  identify  the  aspects  of  the  charge  density  topology  which  most  impact  the  elastic  constants, 
VIP  was  applied  to  the  complete  set  of  parameters  describing  the  CPs  for  the  intermetallic 
compounds.  The  results  of  Figure  8  show  that  all  of  the  CP  descriptors  are  not  equally 
significant  for  the  estimation  of  elastic  constants.  Additionally  different  CP  descriptors  affect  the 
different  constants,  and  therefore  different  topological  features  of  electronic  structure  play  the 
major  role  in  determining  the  mechanical  behaviors.  The  most  significant  descriptors,  which  are 
defined  as  those  having  ‘Variable  Importance’  measures  greater  than  unity,  are  then  used  with 
predictive  informatics  -  partial  least  squares  (PLS)  -  to  develop  the  QSARs.  For  Cn,  17 
variables  meet  the  importance  criteria,  21  variables  for  C12,  and  11  variables  for  C44.  Of  these 
variables,  charge  at  the  cage  critical  points  is  very  significant  for  C44,  shape  at  the  ring  critical 
points  is  very  significant  for  C12,  and  charge  at  the  ring  critical  points  and  the  energy  tenn  related 
to  the  surface  have  significance  for  all  of  the  elastic  constants,  and  therefore  also  anisotropy  and 
Cauchy  pressure.  These  findings  for  the  B2  intermetallics  are  much  different  than  the  descriptor 
(charge  around  bond  critical  point)  previously  discussed  as  controlling  the  FCC  mechanical 
behavior.  While  charge  at  bond  critical  point  is  significant  for  all  three  constants,  it  has 
relatively  minor  effect  as  compared  to  other  descriptors.  A  heat  map  was  created  for  the 
descriptor  data  base  to  further  select  the  relevant  descriptors  based  on  further  selection  criteria 
(Figure  9).  Based  on  this  further  selection  criteria,  all  relevant  descriptors  were  identified  (Table 
5).  We  identify  based  on  the  descriptors  identified  as  detennining  elastic  constants,  Cn  and  C44 
are  detennined  by  the  charge  density  at  the  critical  points,  while  C12  is  detennined  by  the  shape 
at  the  critical  points. 
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Figure  8.  The  importance  measure  of  topological  charge  density  descriptors  with 
respect  to  the  elastic  constant  of  intermetallic  compounds,  which  is  obtained  by 
variable  importance  in  the  projection  (VIP)  scores  for  the  key  predictor  selection. 
The  cut-off  value  of  VIP  is  unity.  1.  //  2.  /?  3.  hbbi  4.  phhi  3.  tan(6)bb 

(=(h_bb  /  kJ  bh)1/2);  6.  A  hb  (=Xi_bb  +  h_bb  +  k3  bb)i  7.  Xi_b;  8.  X2b;  9.  h_bi  10.  pb; 
11.  tan(0)b  (=(h_b  /  ki  b)1/2)i  12.  Aj,  (=X i_b  +  h_b  +  X3j,);  13.  Xi_c;  14.  X2c;  13. 
X3  c;  16.  pc;  17.  A_c  (=X1_c  +  h_c  +  h_c)i  18.  Xi_r;  19.  X2_r;  20.  X3r;  21.  pr;  22.  A_r 
(=h_r  +  h_r  +  hj);  23.  f,  (=pcmin/pbmax);  24.  f 2  (=(pb  -  Pc)2);  25.  f 3  (=(Pb  +  Pc)2f 


Limited  correlation  to  Cij 


Figure  9.  Heat  map  and  dendrograms  for  the  descriptor  space  of  Table  1  and  the 
corresponding  elastic  constants.  This  figure  maps  out  the  correlation  between 
the  descriptors  and  demonstrates  the  complexity  of  the  input  data  base,  thus 
demonstrating  the  necessity  of  data  mining.  From  this  map,  we  identify 
additional  descriptors  that  were  not  identified  through  VIP.  The  additional 
descriptor  potentially  relevant  for  modeling  the  elastic  constants  are  shown  in 
Table  2. 


Table  5.  A  list  of  descriptors  identified  as  most  significant  for  the  respective 
elastic  constants  using  VIP  and  heat  map.  The  descriptors  defined  as  most 
important  are  based  on  their  selection  by  both  methods.  When  considered  in 
magnitude  of  importance,  we  find  that  the  charge  at  the  bond  critical  point  (pf), 
the  shape  at  the  ring  critical  point  0.2  r),  and  the  charge  at  the  cage  critical  point 
(pc)  dictate  Cn,  C12,  and  C44,  respectively. 


Cll 

C12 

C44 

VIP 

Pb,  Pc,  Pr,  fs 

^■2r,  fi,  f3 

Pc,  Pr,  fl,  fi 

Heat 

map 

f3 

Ar,  A.3b 

Pbb,  Pr,  Pb 

fi,  fi,  ^2r, 

Pb,  Pc,  Pr, 

Selected 

Pb,  Pc,  Pr,  f3 

^3b,  Ar 

Pbb,  h,  fi 

(iv)  Informatics  Modelins  of  Critical  Point  Descriptors 

From  the  developed  quantitative  structure-property  relationship  (QSPR),  the  charge  at  the  bond 
critical  points  (pb),  the  shape  at  the  ring  critical  points  (at),  and  the  charge  at  the  cage  critical 
points  (pc)  were  identified  as  the  most  important  features  for  determining  Cll,  C12,  and  C44, 
respectively  (Table  5  and  Figure  10).  While  that  work  accelerated  the  selection  of  chemistries 
for  targeted  elastic  behavior,  it  still  required  the  DFT  calculation  of  the  critical  point  descriptors. 
To  address  this  limitation,  we  use  an  informatics  based  approach  to  “learn”  the  mathematics  of 
the  DFT  calculation  and  predict  these  critical  point  values  for  thousands  of  chemistries  in 
seconds.  Thus  this  work  developed  the  use  of  infonnatics  for  approximating  the  complex 
mathematics  of  DFT,  and  is  used  for  targeted  charge  density  and  elastic  behavior  design.  The 
input  into  the  analysis  is  shown  in  Table  6,  with  the  results  of  the  descriptor  inputs  shown  in 
Figure  1 1 .  From  this  analysis,  we  identified  the  minimum  number  of  elemental  descriptors 
required  to  describe  the  charge  density  topology  and  the  mechanical  response. 
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Figure  10.  The  elastic  constants  (Cll,  C12,  and  C44)  were  modeled  as  a  function  of 
charge  density  critical  points.  Using  a  variable  importance  in  the  projection  approach  to 
identify  the  most  critical  charge  density  descriptors  for  each  elastic  constant,  and  then  a 
partial  least  squares  analysis  of  these  descriptors,  the  controlling  charge  density  features 
were  identified.  For  Cll  charge  at  the  bond  critical  point  was  most  important,  for  Cl 2 
shape  at  the  ring  critical  point  was  most  important,  and  for  C44  charge  at  the  cage 
critical  point  was  most  important.  We  now  predict  these  three  critical  points  as  a 
function  of  the  elemental  components  of  the  intermetal  lies. 


Table  6.  Input  into  the  PLS  analysis.  The  first  three  columns  are  the  DFT  calculated 
critical  point  descriptors. 


Element  A  Element  B 


TMRE 

charge b 

eis2 i 

diarge c 

Pseudopole 
ntial  core 
radii  sum 

R=r  c+r„ 

Pauling 

electronega 

tivity 

Melting 
Point  (K) 

Atomic 

Wei^rt 

Density  @ 
293  K 
(g*m5) 

Pseudopote 

ntial  core 
radii  sum 
R=rt+r„ 

Pauling 

electronega 

tivity 

Melting 
Point  (K) 

Atomic 

Weight 

Density  @ 
293  K 
(g/cm5) 

AgCa 

0.09953 

0,146708 

0.05047 

2.375 

1.93 

1235.08 

107.8682 

10.5 

3 

1 

1112 

40  078 

1.55 

Afll 

0,12461 

0,157876 

0.042866 

2.375 

1.93 

1235  08 

107.8682 

105 

1,585 

2.66 

386.7 

126.90447 

4.93 

AgMo 

0.2S2L39 

0,218766 

0.165314 

2.375 

1.93 

1235.08 

107.8682 

10.5 

2.72 

2.16 

2890 

95.94 

10.22 

Ag  Sc 

0.167966 

0.29196 

0.090283 

2.375 

1.93 

1235.08 

107.8682 

10.5 

2,75 

1.36 

1814 

44.95591 

2.989 

Al(Tl 

0.178147 

0.137461 

0.131953 

2.375 

1  93 

1235.08 

107.8682 

10.5 

2.58 

1.54 

1933 

47.88 

4.54 

AgW 

0.^0246 

0,272557 

0.173413 

2,375 

1.93 

1235.08 

107.8682 

10,5 

2.735 

2.36 

3680 

183.84 

19.3 

AgY 

0.151633 

0.238017 

0,076233 

2.375 

1  93 

1235.08 

107.8682 

10.5 

2.94 

1.22 

1795 

88.90585 

4.469 

KRe 

0.119693 

0.178391 

0.048341 

3,69 

0.82 

336  8 

39  0983 

0.862 

2.68 

1.9 

3453 

186.207 

21.02 

kW 

0.108844 

0.124221 

a  051475 

3.69 

0.82 

336.8 

39.0983 

0,862 

2.735 

2.36 

3680 

183.84 

19.3 

MnMg 

0.231297 

0,074334 

0.034449 

2,22 

1.55 

1517 

54.93805 

7.44 

2.03 

1.31 

922 

24.305 

1.738 

MnSc 

0.213969 

0.191127 

0.139318 

2,22 

1.55 

1517 

54.93805 

7.44 

2.75 

1.36 

1814 

44  95591 

2.989 

MnW 

0. 370212 

0.025284 

0.252933 

2,22 

1.55 

1517 

54.93805 

7.44 

2.735 

236 

3680 

183.84 

19.3 

MnY 

0.191783 

0.21734 

0.114336 

2,22 

1.55 

1517 

54.93805 

7.44 

2.94 

1.22 

1795 

88  90585 

4.469 

Mo  Sc 

0.2)3182 

0,193447 

0.141205 

2.72 

2.16 

2890 

95  94 

10.22 

2.75 

1.36 

1814 

44  95591 

2.989 

MoSr 

0.136147 

0.150869 

0.073498 

2.72 

2.16 

2890 

95  94 

10.22 

3.21 

0.95 

1042 

87.62 

2.54 

MoY 

0.191951 

0.17474 

0.11486 

2.72 

2.16 

2890 

95.94 

10  22 

2,94 

1.22 

1795 

88.90585 

4  469 

NbOs 

0.363461 

0.094869 

0.213102 

2,76 

1.6 

2741 

92.90638 

8.57 

2.65 

2.2 

3327 

190.23 

22.59 

NbRe 

0.360153 

0.012579 

0.215449 

2,76 

1.6 

2741 

92.90638 

8.57 

2,68 

1.9 

3453 

186.207 

21.02 

NbY 

0.170374 

0.123322 

0,114401 

2,76 

1.6 

2741 

92.90638 

8.57 

2.94 

1.22 

1795 

88  90585 

4.469 

NiSc 

0.2)5029 

0,238946 

0.116457 

2.18 

1.91 

1726 

58.6934 

8.902 

2.75 

1.36 

1814 

44  95591 

2.989 

NiY 

0.181377 

0.271383 

0.096722 
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Figure  11.  Weights  of  the  PLS  model,  standardized  based  on  different  units  of  inputs.  As 
can  be  seen,  the  shape  descriptor  is  much  different  than  the  charge  descriptors  (opposite 
sign  for  multiple  descriptors),  while  the  two  charge  descriptors  are  similar,  particularly 
for  the  transition  metal.  To  accelerating  the  modeling  of  the  charge  density  topology,  we 
predict  those  descriptors  with  high  magnitude  weight,  and  remove  those  descriptors  with 
low  weight. 


Using  this  reduced  descriptor  space  from  VIP  analysis  of  the  high-dimensionality  charge 
density  description,  we  applied  partial  least  squares  (PLS),  a  predictive  informatics 
approach  which  accounts  for  multi-collinearity  in  high  dimensional  data,  we  predicted  the 
charge  density  critical  points,  which  we  showed  in  Figure  6  can  then  be  used  to  model  the 
elastic  moduli.  The  results  of  this  prediction  is  shown  in  Figure  12. 
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Figure  12.  Result  of  the  PLS  analysis.  Based  on  the  input  descriptors,  the  charge  at  the 
bond  and  cage  critical  points  were  modeled  with  high  accuracy,  demonstrating  that 
informatics  was  able  to  accurately  model  these  aspects  of  the  charge  density  without 
requiring  any  additional  DFT  calculations.  However,  the  accuracy  for  the  shape 
descriptor  was  insufficiently  accurate,  implying  that  the  controlling  physics  is  different 
than  it  was  for  the  charge. 


In  Figure  12,  we  found  that  the  elemental  descriptors  of  Figure  11  were  sufficient  for 
modeling  the  charge  density  at  the  critical  points.  However,  for  modeling  the  shape  at  the 
critical  points,  these  descriptors  were  insufficient.  Through  further  informatics  analysis, 
we  found  that  to  model  the  shape  at  the  critical  points,  the  elemental  work  functions  must 
be  added.  When  the  work  function  is  included,  the  accuracy  of  the  prediction  of  shape 
improves  by  over  50%  (Figure  13). 
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Figure  13.  Improvement  in  prediction  of  X2 r  by  adding  work  function  into  the  predictor 
set.  Accuracy  improved  from  60%  to  90%.  The  difference  between  charge  density  at  the 
critical  points  and  shape  of  the  charge  density  topology  of  the  critical  point  is  described 
by  the  work  functions  of  the  constituent  elements. 
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We  can  thus  predict  the  elastic  moduli  from  the  critical  point  descriptors  and  we  can  now  model 
the  critical  point  descriptors  from  descriptors  of  the  elements  comprising  the  alloy.  Based  on 
this,  we  can  now  predict  the  electronic  structure  and  mechanical  response  of  ‘virtual’  materials 
based  only  on  information  from  the  periodic  table.  This  work  significantly  accelerates  the 
targeted  design  of  new  materials  based  on  electronic  structure  design. 


(v)Accelerated  Design  Based  on  Density  of  States 

From  the  DOS  spectra,  and  particularly  from  the  symmetry  induced  d-orbital  splitting,  of  each 
transition  metal  in  each  of  bcc,  fee,  and  hep  structure  (provided  by  the  CSM  group),  we  have 
applied  informatics  methods  to  these  d-orbital  splittings  to  predict  the  ground  state  structure  of 
the  transition  metals  (Figure  14).  The  prediction  based  on  pattern  recognition  and  partial  least 
squares  regression  is  found  to  work  very  well  for  this  class  of  material. 
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Figure  14. Informatics  prediction  of  ground  state  ciystal  structure  of  second  row 
transition  metals.  The  model  predicts  crystal  structure  from  a  DOS  input  of  each 
element  in  three  structures:  fee,  bcc,  and  hep  structures.  The  red  squares  are  elements 
used  in  training  the  model,  while  the  blue  diamonds  are  elements  which  were  not 
included  in  building  the  model  but  were  used  instead  for  external  validation. 


Building  on  this  structural  modeling  from  DOS  curves,  we  use  the  approach  to  model 
unknown  DOS  spectra  using  artificial  neural  networks  (ANN).  Using  an  input  of  DOS  of 
Fe,  Ni,  Pd,  Fe3Ni,  FeNfi,  Fe3Pd,  and  FePd3,  a  model  linking  these  elemental  DOS  to 
multi-component  alloys  was  developed.  This  model  is  then  tested  using  other  known 
alloy  DOS  not  used  in  model  development.  The  necessary  input  to  model  these  systems 
are  the  DOS  spectra  for  Fe  (FCC),  Pd  (BCC)  and  Ni  (BCC).  Also  DOS  spectra  of  two 
alloys  in  each  system  were  introduced  so  that  the  learning  technique  can  capture  the 
nature  of  the  interactions  of  elements  within  alloys.  Along  with  the  DOS  spectra  of  the 
elements  of  the  binary  systems,  the  DOS  spectra  of  two  chemistries  of  each  binary  system 
were  used  to  train  the  model.  The  model  was  tested  for  FePd  and  FeNi  (Figure  15),  which 
were  not  included  in  any  form  in  training  of  the  model,  but  were  used  solely  for  testing 
the  accuracy  of  the  model.  Figure  15  shows  these  DOS  spectra  with  the  prediction 
compared  to  DFT.  The  model  for  predicting  DOS  of  alloy  systems  shows  a  level  of 
accuracy,  with  the  general  shape  of  the  DOS  curve  captured  but  the  fine  structure  is  not 
fully  captured.  The  model  therefore  links  the  elemental  DOS  with  the  alloy  DOS.  Of 
particular  note  is  that  crystal  structure  is  not  explicitly  defined  in  this  analysis.  The  only 
description  of  the  alloy  necessary  is  the  stoichiometry.  For  the  results  shown,  the  input 
for  FexNiy  was  x  Fe  and  y  Ni,  while  similarly  for  FexPdy  the  input  was  x  Fe  and  y  Pd  . 
However,  this  prediction  quality  is  comparable  with  the  other  techniques  for  modeling 
multi-component  systems  without  need  of  describing  the  atomic  interactions,  but  with 
more  applicability  by  not  requiring  an  input  Hamiltonian,  and  is  sufficient  for 
representing  the  significant  features  of  the  DOS  curve.  Additionally,  this  level  of 
accuracy  is  particularly  notable  given  the  limited  amount  of  input  information  needed.  As 


most  elemental  DOS  are  known,  this  approach  then  allows  us  to  model  DOS  of  any 
chemistry  and  stoichiometry. 

In  these  results,  the  accuracy  of  the  model  is  comparable  to  that  for  the  training 
systems,  indicating  that  our  model  is  not  over-fitting  the  data.  The  crystal  structures  for 
these  systems  were  not  input  into  the  model,  but  we  still  capture  the  interactions 
accurately.  For  FeNi  and  FePd,  the  crystal  structure  type  is  Ll0,  while  for  FeX3  and  Fe^X 
(X=Pd,  Ni),  the  crystal  structure  type  is  LI2.  Capturing  the  general  shape  of  the  DOS 
curves  indicates  that  ANN  is  able  to  extract  the  electronic  interactions  from  the  elemental 
DOS,  and  is  able  to  represent  the  form  of  these  electronic  interactions  in  multi-component 
systems. 


Figure  15.  Testing  of  DOS  modeled  for  FePd  and  FeNi  using  the  trained  ANN.  These  two 
chemistries  were  not  included  in  the  development  of  the  ANN  model,  and  serve  as 
validation  sets  of  the  approach.  The  accuracy  of  the  test  DOS  is  comparable  with  the 
training  DOS,  indicating  that  we  have  avoided  an  over-fitting  of  the  data,  and  the  model 
is  applicable  for  other  stoichiometries  of  the  elements. 


Using  the  developed  logic,  the  bulk  modulus  of  the  training  and  test  systems  were 
predicted  (Figure  16).  From  this  figure,  we  find  that  the  logic  is  capable  of  predicting 
bulk  modulus  from  the  DOS  spectra  with  high  accuracy.  Additionally,  the  model  predicts 
bulk  modulus  well  for  both  the  training  and  test  systems,  indicating  that  our  model  is  not 
over-fitting  the  training  data.  The  input  for  the  test  systems  was  only  ANN  modeled 


DOS,  with  no  input  of  bulk  modulus.  The  accuracy  of  this  result  validates  that  the  level 
of  detail  of  the  ANN  model  is  sufficient  for  describing  the  alloys.  This  work,  combined 
with  the  work  on  charge  density  topology,  shows  the  value  of  informatics  analysis  of  ab 
initio  calculated  data. 
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Figure  16.  The  result  of  the  approach  for  predicting  bidk  modulus  from  an  input  of  alloy 
DOS  spectra.  The  diamonds  correspond  to  DOS  calculated  by  DFT  while  the  squares 
correspond  to  DOS  modeled  via  ANN.  The  bulk  modulus  of  the  training  systems  was  used 
to  train  the  model,  while  the  bulk  modulus  of  the  test  systems  was  not  used  in  any  form  in 
developing  the  model.  The  accuracy  of  the  results  demonstrates  that  bulk  modulus  is 
clearly  represented  within  the  DOS  spectra,  and  that  it  can  be  quantitatively  extracted  via 
statistical  learning.  Additionally,  the  accuracy  of  the  bulk  modulus  predicted  for  the  test 
systems  indicates  that  the  ANN  approach  for  modeling  DOS  is  sufficiently  detailed  and 
accurate. 
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